############## 
############## 
############## Models
remove(list = ls())

base::library(texreg)
base::library(conflicted)
base::library(tidyverse)
conflict_prefer("filter","dplyr")
base::library(ggplot2)
base::library(dplyr)
base::library(here)
conflict_prefer("here", "here")
base::library(systemfit)
base::library(stargazer)
base::library(DataCombine)
base::library(fixest)
base::library(ggeffects)

load(here("Data", "Model_Data.RData"))
glimpse(final_data)

summary(final_data$death_rate)
table(final_data$year)

final_data <- final_data %>% 
  mutate(time_trend = ifelse(year == 2004, 1, NA),
         time_trend = ifelse(year == 2008, 2, time_trend),
         time_trend = ifelse(year == 2012, 3, time_trend),
         time_trend = ifelse(year == 2016, 4, time_trend),
         time_trend = ifelse(year == 2020, 5, time_trend)) 

glimpse(final_data)


final_data <- final_data %>% 
  mutate(quartile = ifelse(death_rate <= 7.351, 1, NA),
         quartile = ifelse(death_rate > 7.351, 2, quartile),
         quartile = ifelse(death_rate > 11.522, 3, quartile),
         quartile = ifelse(death_rate > 18.295, 4, quartile)) %>% 
  mutate(quartile = factor(quartile, levels = c(1, 2, 3, 4)))
table(final_data$quartile)


eq1_d <- dem_rep ~ quartile + unemp_rate + 
  log(defl_pcincome) + 
  share_black + share_hisp + share_asian + share_young + share_elder + 
  log1p(share_college) + log1p(share_manuf) +
  log(tot_pop) + as.factor(rural_urban) + state_partiesdiff +
  state + rep_inc + incumbent_cand + lag_dem_rep + time_trend
eq2_d <- abs_rep ~ quartile + unemp_rate + 
  log(defl_pcincome) + 
  share_black + share_hisp + share_asian + share_young + share_elder + 
  log1p(share_college) + log1p(share_manuf) +
  log(tot_pop) + as.factor(rural_urban) + state_partiesdiff +
  state + rep_inc + incumbent_cand + lag_abs_rep + time_trend
Syst_death <- list(eq1_d, eq2_d)
model_d <- systemfit(Syst_death, method = "SUR", data = final_data)
summary(model_d)


summary(ols_model1 <- lm(rep_share2 ~ quartile + unemp_rate + 
                             log(defl_pcincome) + 
                             share_black + share_hisp + share_asian + share_young + share_elder + 
                             log1p(share_college) + log1p(share_manuf) +
                             log(tot_pop) + as.factor(rural_urban) +
                             state + state_partiesdiff + rep_inc + incumbent_cand + lag_rep_share2 + time_trend,
                           data = final_data))

summary(ols_model2 <- lm(rep_share2 ~ death_rate + I(death_rate^2) + unemp_rate + 
                           log(defl_pcincome) + 
                           share_black + share_hisp + share_asian + share_young + share_elder + 
                           log1p(share_college) + log1p(share_manuf) +
                           log(tot_pop) + as.factor(rural_urban) +
                           state + state_partiesdiff + rep_inc + incumbent_cand + lag_rep_share2 + time_trend,
                         data = final_data))

summary(ols_model3 <- lm(dem_rep ~ quartile + unemp_rate +
                           log(defl_pcincome) + 
                           share_black + share_hisp + share_asian + share_young + share_elder + 
                           log1p(share_college) + log1p(share_manuf) +
                           log(tot_pop) + as.factor(rural_urban) +
                           state + state_partiesdiff + rep_inc + incumbent_cand + lag_dem_rep + time_trend,
                         data = final_data))

summary(ols_model4 <- lm(abs_rep ~ quartile +unemp_rate +
                           log(defl_pcincome) + 
                           share_black + share_hisp + share_asian + share_young + share_elder + 
                           log1p(share_college) + log1p(share_manuf) +
                           log(tot_pop) + as.factor(rural_urban) +
                           state + state_partiesdiff + rep_inc + incumbent_cand + lag_abs_rep + time_trend,
                         data = final_data))


stargazer(ols_model1, ols_model2, ols_model3, ols_model4,
          type = "latex",
          title = "OLS and Fixed Effect Estimates",
          dep.var.labels = c("Republican Vote Share", "Republican Vote Share", "log(Dem/Rep)", "log(Abs/Rep)"), 
          no.space = FALSE, single.row = T,
          covariate.labels = c("Overdose Death Rate (2nd Quartile)",
                               "Overdose Death Rate (3rd Quartile)",
                               "Overdose Death Rate (4th Quartile)",
                               "Overdose Death Rate (level)",
                               "Overdose Death Rate (cube)",
                               "Local Unemployment Rate",
                               "Log PC Income",
                               "Share Black",
                               "Share Hispanic",
                               "Share Asian",
                               "Share Young Voters",
                               "Share Elder Voters",
                               "Log Share College Degree",
                               "Log Share Manufacture",
                               "Log Population",
                               "Close Election (State)",
                               "Republican Presidency",
                               "Incumbent Candidate",
                               "Lagged Republican Vote Share"),
          omit = c("rural_urban", "state"),
          omit.labels = c("Rural-Urban Code",
                          "State Fixed Effects"))


model1 <- ggpredict(ols_model2, terms = "death_rate", ci.lvl = 0.90)
glimpse(model1)
ggplot(model1, aes(x = x, y = predicted)) +
  geom_line(linewidth = .7, colour = "midnightblue") + 
  theme_bw() +
  geom_ribbon(aes(ymin = conf.low, ymax = conf.high), alpha = .1,
              fill = "#78b3c0") +
  scale_x_continuous(breaks = c(20, 40, 60, 80, 100), lim = c(0, 101)) +
  ggtitle(NULL) +
  xlab("Death Rates") +
  ylab("Predicted Republican Vote Share") +
  theme(axis.text.x = element_text(size = 13#, 
                                   #hjust = 1, 
                                   #angle = 15
  ),
  axis.text.y = element_text(size = 12),
  axis.title.y = element_text(size = 13),
  legend.text = element_text(size = 13),
  legend.title = element_text(size = 12),
  legend.position = "right",
  legend.direction = "vertical",
  title = element_text(size = 14))



